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ABSTRACT 

We investigate the production of cosmic ray (CR) protons at cosmological shocks by performing, for the 
first time, numerical simulations of large scale structure formation that include directly the acceleration, 
transport and energy losses of the high energy particles. CRs are injected at shocks according to the 
thermal leakage model and, thereafter, accelerated to a power-law distribution as indicated by the test 
particle limit of the diffusive shock acceleration theory. The evolution of the CR protons accounts for 
losses due to adiabatic expansion/compression, Coulomb collisions and inelastic p-p scattering. Our 
results suggest that CR protons produced at shocks formed in association with the process of large 
scale structure formation could amount to a substantial fraction of the total pressure in the intra-cluster 
medium. Their presence should be easily revealed by GLAST through detection of 7-ray flux from the 
decay of ir° produced in inelastic p-p collisions of such CR protons with nuclei of the intra-cluster gas. 
This measurement will allow a direct determination of the CR pressure contribution in the intra-cluster 
medium. We also find that the spatial distribution of CR is typically more irregular than that of the 
thermal gas because it is more influenced by the underlying distribution of shocks. This feature is 
reflected in the appearance of our 7-ray synthetic images. Finally, the average CR pressure distribution 
appears statistically slightly more extended than the thermal pressure. 

Subject headings: acceleration of particles — cosmology: large-scale structure of universe — gamma 
rays: theory — methods: numerical — shock waves — X-rays: galaxies: clusters 



1. INTRODUCTION 

Clusters of galaxies are the largest bound objects in the 
universe and prove invaluable for investigations of cosmo- 
logical interests. The statistics of cluster masses and their 
dynamical properties, including, for instance, the relative 
proportions of baryonic and non-baryonic matter, are com- 
monly used to test basic cosmological models (e.g., Bah- 
call 1999, and references therein). While galaxies are the 
most obvious constituents of clusters in visible light, most 
of the cluster mass is non-baryonic, and even the bary- 
onic matter is primarily contained within the diffuse intra- 
cluster medium (ICM), rather than in galaxies. The tem- 
perature and density distribution of the ICM gas directly 
reflect the dynamical state of clusters, a topic that has re- 
ceived much attention recently. While the ICM of clusters 
sometimes appears relaxed, it is often the case that high 
speed flows are present, demonstrating that cluster envi- 
ronments can be violent (e.g., Markevitch et al. 1999, and 
references therein). 

The likely existence of strong "accretion" shocks sev- 
eral Mpc's from cluster cores developed in the course of 
large-scale structure formation has been recognized for a 
long time (Sunyaev & Zel'dovich 1972; Bertschinger 1985; 
Ryu & Kang 1997; Quilis et al. 1988). Such shocks are 
responsible for the heating of the ICM up to temperatures 
of order 10 — 10 8 K. However, cosmic structure formation 



simulations have demonstrated, in addition to accretion 
shocks and discrete merger shocks, the existence of some- 
what weaker shocks "internal" to the ICM that are very 
common and complex (Miniati et al. 2000). Since clus- 
ters tend to form at the intersections of cosmic filaments, 
they accrete matter in unsteady and non-isotropic patterns 
through large scale flows propagating down filaments and 
producing shocks as they impact the ICM. When cluster 
mergers take place, the accretion shocks associated with 
the individual clusters add to the shocks that form in di- 
rect response to the merger. The net result of all of this is 
a rich web of relatively weak shocks, which often penetrate 
into the inner regions of the clusters (Miniati et al. 2000). 
Shocks resulting from discrete cluster merger events have 
already been identified by the observation of temperature 
structure in the ICM (e.g., Markevitch et al. 1999, and 
references therein). Such shocks have also been claimed as 
the acceleration sites for relativistic electrons responsible 
for the non-thermal emission observed from clusters in the 
radio, hard X-ray and extreme ultra-violet (e.g., Takizawa 
& Naito 2000; Enfilin et al. 1998; Roettiger et al. 1999, and 
references therein). 

Magnetic fields are commonly observed in the large scale 
structures (e.g., Kronberg 1994). They may have been 
seeded at shocks in the course of structure formation and 
amplified up to /iG level in clusters and, perhaps, also 
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in filaments and super-clusters (Kulsrud et al. 1997; Ryu 
et al. 1998). Because shock waves in the presence of even 
modest magnetic fields are sites of efficient cosmic ray 
(CR) acceleration (e.g., Drury 1983), structure formation 
might imply copious generation of high energy particles, 
including both protons and electrons. In fact, according to 
diffusive shock acceleration theory (Drury 1983), as much 
as several tens of percent of the kinetic energy of the bulk 
flow associated with the shock can be converted into CR 
protons (Eichler 1979; Berezhko & Ellison 1999). Given 
their huge size and long durability, large scale structure 
shocks have also been suggested as possible sources of the 
very high energy CR protons up to a few x 10 19 eV (Kang 
et al. 1996, 1997). 

The overall energetics of the "cosmic" (accretion and 
internal) shocks is generally consistent with the produc- 
tion of CRs containing a significant energy fraction. Typ- 
ical flow speeds in and around clusters will be ~ 
(2G Md/Rci) 1 / 2 ~ 2 x 10 3 km s"\ leading to an avail- 
able power for CR production at accretion shocks ~ 
p b v 3 f R 2 cl ~ 10 46 erg s" 1 using M cl ~ 1O 15 M and R cl ~ 2 
Mpc. According to Miniati et al. (2000), where the statis- 
tics of cosmic shocks in both SCDM and ACDM cosmolo- 
gies are explored, accretion shocks appear to be less impor- 
tant as potential sources of CRs than internal shocks, de- 
spite the typically greater strength of the accretion shocks. 
The reason is that internal shocks repeatedly process the 
ICM material, whereas the accretion shocks do it only once 
with low density background material (cf. Miniati et al. 
2000, for $e in clusters estimated from simulations). 

Observations of radio emission from CR electrons as 
well as radiation excess in the hard X-ray and possibly 
EUV bands, have recently stimulated much discussion 
about cluster physics (e.g., Sarazin 1999, and references 
therein). They have provided information regarding the 
energy density of CR electrons. CR protons produce 7- 
rays through n° decay following inelastic collisions with 
gas nuclei. While such 7-rays have not yet been detected 
from clusters (Sreekumar et al. 1996), recent estimates 
have shown that 7-ray fluxes from the nearest rich clus- 
ters, such as Coma, are within the range of what may 
be detected by the next generation of 7-ray observatories 
(Colafrancesco & Blasi 1998; Blasi 1999; Dolag & Enfilin 
2000). Their detection will provide essential information 
about the presence of and amount of energy carried by CR 
protons in the ICM (cf. §3.3). 

Rclativistic protons below the "GZK" energy threshold 
for photo-pion production due to interaction with the cos- 
mic microwave background photons (i.e., E < 10 9 ' 5 GeV) 
do not suffer significant energy losses in cluster environ- 
ments during a Hubble time (Berezinsky et al. 1997). In 
addition, up to somewhat lower energies (~ 10 6 GeV), even 
conservative estimates of diffusion rates would confine CRs 
within clusters since their formation (Volk et al. 1996; 
Berezinsky et al. 1997). Therefore, CR protons, once intro- 
duced, should accumulate in clusters, with the possibility 
of impacting on a wide range of issues. Some topics that 
could be impacted include cluster formation and evolution, 
as well as cluster mass estimates based on the assumption 
of ICM Hydrostatic equilibrium. The dynamics of cooling 
flows also would obviously be affected. 

The analysis of Miniati et al. (2000) showed that the 



most common shocks in the ICM have typical Mach num- 
bers less than 10, with a peak around M ~ 4 — 5. That 
is significant, because such shocks are strong enough to 
transfer as much as 20 — 30% of the bulk kinetic energy 
into CR pressure, but are only mildly modified by the CR 
back-reaction (Jones et al. 2000). Thus, the test-particle 
approximation, in which such dynamical feed-back is ig- 
nored, should be a reasonable, physically justified assump- 
tion to begin investigating the production of CR at cosmic 
shocks. 

In this paper we investigate the acceleration of CR pro- 
tons at cosmological shocks by means of numerical calcu- 
lations. For the first time the CR population is directly 
included in the computation with particle injection, ac- 
celeration and energy losses calculated in accord with the 
properties of the local enviroment in which the particles 
are propagating. Here, our focus is CR protons, while CR 
electrons will be discussed in a companion paper (Miniati 
et al. 2001). There are additional sources of CRs in clus- 
ters, of course, such as active galaxies (Enfilin et al. 1997; 
Berezinsky et al. 1997). We do not attempt to include 
them in our current simulations, since our goal is to un- 
derstand the role of structure shocks. However, we do call 
attention in our discussions to some expected differences 
between shock CR sources and point sources as appropri- 
ate. The results of our modeling efforts should provide 
some initial clue as to how these different sources can be 
distinguished observationally. 

The paper is organized as follows. In §2 we outline the 
numerical methods adopted for our study. In §3 the re- 
sults are presented. A discussion on the implication of 
the results of this paper is given in §4, whereas the main 
conclusions are summarized in §5. 

2. NUMERICAL SIMULATIONS 

2.1. Cosmological Hydrodynamic Simulations 

For the numerical calculations we employed an Eulerian 
"TVD" hydro + N-body cosmological code (Ryu et al. 
1993). Since the computation involves a new quantity 
never simulated before in this context, i.e., CRs, we have 
decided to begin the study from the simpler case of the 
standard cold dark matter (SCDM) model, leaving the 
currently more favored CDM + A model as the natural 
follow-up step for future work. Although it is well known 
that SCDM is not a viable model anymore (e.g., Ostrikcr 
1993), we have chosen the key cosmological parameters 
so that properties of the simulated collapsed objects are 
consistent with observations, thus allowing assessments of 
their general characteristics. For instance, we adopted rms 
density fluctuations on a scale of 8/i _1 Mpc to be defined 
by os = 0.6, which is incompatible with COBE results 
and the SCDM universe, yet induces the emergence of a 
reasonable population of collapsed objects in simulations 
of large scale structure formation (Ostriker & Cen 1996). 
We have also adopted the following key parameters: spec- 
tral index for the initial power spectrum of perturbations 
n = 1, normalized Hubble constant h = H o /(100 km s _1 
Mpc" 1 ) = 0.5, total mass density Qm — 1, and bary- 
onic fraction £1^ = 0.13. In addition, we use a standard 
metal composition with hydrogen and Helium mass frac- 
tions fn = 0.76 and fn e = 0.24 respectively. Thus, for 
a fully ionized gas the mean molecular weight used in the 
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temperature definition is p, — 0.59. 

In order to simulate a cosmological volume large enough 
to contain groups/clusters with a sufficient resolution, we 
select a cubic comoving region of size 50 /i _1 Mpc and 
use 256 3 cells for baryonic matter and 128 3 dark matter 
(DM) particles. This corresponds to a spatial resolution 
of ^200/i _1 kpc. A few comments about the effects due to 
finite numerical resolution are appropriate here in order to 
define the scope of our findings. In general a coarse grid 
limits the structures that can form during the evolution 
of the simulated systems. This implies first that density 
peaks are smoothed out while masses of the structures are 
conserved. As a consequence, quantities such as the X-ray 
and 7-ray luminosity, which depend on the square of the 
density, will be reduced. The effect is stronger for lower 
temperature groups/clusters which have similar structures 
as the larger clusters but smaller physical scales. Thus, be- 
cause of resolution effects, these types of emissions (X-ray, 
7-ray) will be systematically underestimated and will lead 
to steeper intra-cluster temperature dependences in our 
numerical calculation. Previous numerical studies carried 
out to test the performance of the hydrodynamic part of 
the code employed here indicate that the X-ray thermal 
emission (oc n 2 gas ) is underestimated by a factor of a few 
(Cen & Ostrikcr 1999). Secondly, although shocks are cap- 
tured cleanly within only a few computational zones, that 
still amounts to a fair fraction of the cluster size. This in- 
troduces an uncertainty in the location of the shocks and 
reduces both the shock surface extension and complexity. 
However, the total flux of kinetic energy through shocks 
should not be affected, as indirectly attested by the fact 
that the computed intra-cluster temperatures are quite ac- 
curate (Kang et al. 1994; Frenk et al. 1999). Thus we ex- 
pect our results to be physically correct, although further 
work is required in order to achieve high quantitative ac- 
curacy. Since this is the first attempt to investigate such 
a problem, the level of accuracy characterizing our simu- 
lation should be sufficient enough to explore qualitatively 
the physical impact that CRs may have in cosmological 
environment and to provide a preliminary assessment of 
their observability. 

2.2. Cosmic Ray Injection and Acceleration 

The evolution of the CR population in the simulation is 
computed via passive quantities by the code COSMOCR 
(Miniati 2001). In the following sections we provide a brief 
description of the physical processes included in this code, 
i.e., CR injection at shocks and spatial transport and en- 
ergy losses. 

In the calculation CRs are injected at shocks according 
to the "thermal leakage" model (e.g., Ellison & Eichler 
1984; Kang & Jones 1995). In this model the post-shock 
gas is assumed to have mostly thermalized to a Maxwellian 
distribution, /(p)Maxweii, characterized by the downstream 
temperature, T s hock- Thermal protons in the high energy 
tail of such a Maxwellian distribution can escape back up- 
stream of the shock if their speeds are sufficient to al- 
low them to avoid being trapped by the plasma waves 
that moderate the shock (Malkov & Volk 1995). Those 
protons are injected into the diffusive shock acceleration 
mechanism and can be accelerated to high energies. In 
the present calculation, the momentum threshold for in- 



jection, Pi n j, is set to a few times the peak thermal value, 
i-e-, 

Pinj = ci 2^/m p fc B r s hock (2-1) 
where m p is the proton mass, ks is the Boltzmann's con- 
stant, T s hock is the post-shock gas temperature, and Ci is 
a parameter which regulates the number of injected par- 
ticles (see below). This limit was chosen to be consistent 
with more detailed, nonlinear CR acceleration simulation 
results described at the end of this subsection. In the test- 
particle limit adopted here, the diffusively accelerated CRs 
emerging from a shock are characterized by a power law 
distribution function given by 

( p Y q 

/(?>)shock = /(Pinj)Maxwel] I ) • (2-2) 

\Pinj J 

extending from pi n j to p m ax- Here, the log-slope is deter- 
mined by the shock strength, i.e., q = 3r/(r — 1) (where 
r is the shock compression ratio); and the normaliza- 
tion is given by the value of the Maxwellian gas distri- 
bution at momentum Pi n j- Thus the thermal distribution, 
/(p)Maxweii, and CR distributions, /(p) s h oc k, join smoothly 
in terms of the momentum coordinate. This power-law CR 
distribution is assigned to each grid cell that is identified 
as "being shocked" within a time step in the numerical 
simulation. The physical upper bound to the CR momen- 
tum distribution is determined by several factors, includ- 
ing the time available for acceleration compared to the 
mean time for particles to re-cross the shock due to the 
competition between wave scattering and advection, the 
extent of a shock compared to particle scattering lengths 
and any competition from energy losses during accelera- 
tion. For parameters appropriate to groups/clusters we 
expect the acceleration to proceed relatively quickly up 
to momenta at least as great as 10 6 GeV/c. CRs with 
even higher energy can be produced in principle (Kang 
et al. 1996). Conservative estimates indicate that these 
very high energy CRs can diffuse out of clusters carrying 
away some energy. That would affect our results only if the 
spectra of the accelerated CRs are significantly flattened 
with respect to the test particle limit above our adopted 
momentum upper limit. However, this type of behavior 
typically is expected only for CR dominated and strongly 
modified shocks. From the observed properties of the in- 
tercluster medium, where most of the pressure is thermal, 
most likely that type of shock does not occur. 

In the simulation we assume the power law is formed 
within one dynamical time step up to p m ax — 10 6 GeV/c, 
so that spatial diffusion of CRs can be neglected and the 
computational cost much reduced. To follow the evolu- 
tion of the CR distribution in detail from injection to 
Pmax would be completely impractical, since it would ne- 
cessitate numerical resolution on the scale of the physical 
thickness of the shocks (Jones et al. 1999). Similarly, since 
spatial diffusion of CRs away from shocks is likely to be 
slow below p max , it is neglected there, but COSMOCR 
does include adiabatic energy changes, as well as the en- 
ergy losses from Coulomb and inelastic p-p collisions with 
the thermal ICM. To do this a Fokkcr-Planck equation is 
solved that has been integrated over finite momentum bins 
to take advantage of the near-power-law form of the CR 
momentum distribution, f(p). In effect, the momentum 
space is divided into 8 logarithmically equidistant inter- 
vals, bounded by pi,...p$,pg = p max , which we refer to 
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here as momentum bins. Within each momentum bin, j, 
we assume /(x i; p) oc p~u( Xi ) , where qj(x.i)) is determined 
self-consistently from n(xi,pj) defined below and the re- 
quired continuity of f(p). At each computational spatial 
grid point, Xj, and for each momentum bin, j, we define 
the number density as 

47T [ Pj + 1 

n{xi,Pj) = it / f(*i,p)p dp. (2-3) 

For a full description of the code COSMOCR we refer to 
Miniati (2001, but see also Jones et al. 1999). 

Before concluding this section, we return for a moment 
to the "injection" parameter c\ which deserves some fur- 
ther comments. As already pointed out the value of c\ 
determines the fraction rji n j of particles in the post-shock 
gas with density n<i that are injected into the CRs as fol- 
lows (Miniati 2001): 

where, obviously, n in j is the number of injected particles. 
In addition, we find the ratio of the CR pressure to ram 
pressure {p\U 2 s , which supplies the energy for the cosmic 
rays) to be (Miniati 2001) 



2 / m p c 

\ Pinj 



3-q 



4-g 



- 1 



A-q 



(2-5) 

where c and u s are the light and shock speed, respectively 
and where we neglected the non-relativistic contribution 
of the CR pressure (this is justified since q ~ 4 so that 
most of the CR pressure is produced by relativistic parti- 
cles; see §3). Both observational and theoretical studies of 
diffusive shock acceleration suggest that canonical values 
of ci should be around 2.3 — 2.5, corresponding to the value 
of rji n j ranging between a few xlO -3 to 10~ 4 (Lee 1982; 
Quest 1988; Kang & Jones 1995). According to the sim- 
ulations by Gieseler et al. (2000) where a self-consistent 
injection treatment based on the plasma physical model 
of Malkov & Volk (1995) is adopted, the above injection 
parameter is in fact very reasonable. For the flow param- 
eters relevant for the cosmic shocks in our simulations, we 
find that a value of c\ = 2.6 produces an injection effi- 
ciency rjinj and a post-shock CR pressure consistent with 
the saturation value obtained from numerical studies of 
shock acceleration in which the back reaction of the parti- 
cles is accounted for (Berezhko et al. 1995). Note that, for 
this reason, our value of c\ is somewhat larger than the 
canonical values due to the test-particle treatment. For a 
shock with Mach number M <~ 4 and speed u s <~ 10 3 km 
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theory (Jones et al. 2000). So our choice of c\ is also con- 
sistent with our assumption that the CR acceleration can 
be treated by the test-particle theory. 

2.3. Extracting Global Properties of Groups/ Clusters 

After the calculation was completed, the simulated 
groups/clusters have been identified by the DM-based 
"spherical over-density" method described in (Lacey & 
Cole 1994). The details of the group/cluster identification 
procedure can be found in Miniati et al. (2000). Global 
group/cluster properties, such as core temperature, aver- 
age pressure, emissivity at various wavelengths etc., were 
calculated by averaging or integrating the quantities over 
the group/cluster volume. These global properties have 
then been studied by means of correlation plots in order 
to make predictions about the quantities under investiga- 
tion (and often yet to be measured) in terms of the well 
established ones (see §3.1-3.3). We point out from the 
outset that, because of the relatively small computational 
box, the temperature of the simulated collapsed objects 
only ranges between 0.3 and 3 keV. Nevertheless, after 
determining the temperature dependence of the various 
properties of interest we extrapolate their values beyond 
these temperature limits and make estimates even for rich 
clusters such as Coma which are easier to observe. So long 
as there are no important scales involved, these extrapo- 
lations should be reliable. 

Two-dimensional projections of individual group/cluster 
structures have also been constructed from the data-set, 
either as slices through the simulated volume or as syn- 
thetic images of cluster emissions. The synthetic images, 
computed here in the X-ray and 7-ray bands, are pro- 
duced by means of a projection code (Trcgillis 2001) that 
integrates the "emissivity" along the line of sight in the 
optically thin plasma approximation. These images al- 
low a more in-depth inspection of the spatial distribution 
of the quantities of interests, but for space reasons are 
limited here to only a few examples (see §3.4). 

In general, the 7-ray flux and the surface brightness for 
the synthetic images have been calculated by arbitrar- 
ily setting the groups/clusters to a luminosity distance 
of about 70 /i _1 Mpc (i.e., z — 0.023) corresponding to 
the Coma cluster (Abell 1656). Since our grid resolution 
amounts to ~ 200 /i _1 kpc, at this distance the minimum 
size of a pixel of the synthetic image corresponds to 9.8 
arc-min square. 

3. RESULTS 

3.1. Cosmic Ray Energy Content 
The CR pressure at Xj is defined by 



P\u% \ 4 J \10 3 km s~ 

(2-6) 

Since most of the flow kinetic energy is processed by cosmic 
shocks with M <~ 4 — 5, we expect from experience with 
detailed CR shock simulations that up to 10-30 % of it will 
be converted into CR energy with the above value of c\. 
Detailed simulations also show that modifications to such 
shocks are small enough that the form of the CR spec- 
trum is not substantially changed from the test-particle 
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with f(xi,p) reconstructed from n(xi,p) and q(xi,p) as de- 
scribed in §2.2. From our simulations we find that n(xi,p) 
has a strong spatial dependence, whereas q(xi,p) assumes 
a relatively narrow range of values, mostly between 4.01 
and 4.2. The thermal pressure obeys by the equation of 
state for an ideal gas 



Pth(xi) = n tot (xi) k B T(xi), 



(3-2) 
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where n to t = n lon + n e = (2f H + 3f H e/4:)p g as/m p is the 
gas number density inclusive of both ions and electrons, 
and T the gas temperature. From the thermal and CR 
pressures defined at each cell, we calculate the mean ther- 
mal and CR pressure of groups / clusters within a sphere of 
radius R c i ~ 0.5 /i _1 Mpc from the cluster center as 

where the summation over i extends to the groups/clusters 
volume V = and the weight function W{ is given 

by the portion of each computational cell within V. Given 
our resolution, the volume within a radius of 0.5 h~ 1 Mpc 
typically includes about 65 computational cells. We also 
calculate the groups/clusters core temperature, T x as a 
volume averaged ICM temperature within the same vol- 
ume. 

The first important result of this study is illustrated 
in Fig. 1. There we plot the ratio (P cr /Pth)d for 
groups/clusters at the current epoch (i.e., z = 0) as a 
function of the core temperature, T x . The values of this 
ratio, (P cr /Pth)ci ^ {E cr /2E th ) c i, where E stands for en- 
ergy density, offers a first-order indication of the relative 
importance of the two components for the dynamics of 
groups/clusters. From Fig. 1, we can read that a signif- 
icant fraction (up to ~ 45%) of the total pressure inside 
today's groups/clusters could be borne by CRs. We note 
here that the actual content of CR pressure (and energy 
density) depends on the injection parameter for which we 
have made an educated estimate based on published stud- 
ies related to the theory of diffusive shock acceleration. 
As already pointed out in §§1 and 2.2 here we are only al- 
lowed a simplified treatment of the injection mechanism 
at shocks. Therefore, the result in Fig. 1 should not 
be interpreted as a precise estimate of the CR content 
inside groups/clusters of galaxies. Rather, it provides a 
qualitative, yet sound, insight that CRs might be quite 
important for the dynamics of those objects. Consider- 
ing the difficulties in following the physics of CR accel- 
eration self-consistently in multi-dimensional simulations, 
the quantitative estimate of the total CR content needs 
to be done through measurements of 7-ray fluxes from 
groups/clusters, as we shall see below. 

The ratio of pressures plotted in Fig. 1 does not show 
any particular trend except for a slight reduction toward 
higher temperatures. This might be due to our injection 
model which, according to eq. 2-5, produces a higher ra- 
tio of CR to thermal pressure for shocks with smaller ve- 
locities (u s ) and similar Mach numbers (yielding similar 
q(M)); i.e., cooler preshock gas. It is possible that such 
shocks occur in cooler groups/clusters, characterized by 
lower accretion velocities and similar pre-shock tempera- 
tures. However, the trend in Fig. 1 could also be due to 
adiabatic compression inside the cluster, which increases 
the thermal pressure at a higher rate than the CR one. At 
least we can be sure that the apparent scatter there is in 
part a reflection of the diverse CR acceleration histories of 
groups/clusters. In addition, part of this scatter can also 
be due to the different spatial distribution of the thermal 
and CR components as we shall see in §3.2 and 3.4. 

3.2. Spatial Distribution of Thermal and Cosmic Ray 

Pressure 



Another feature of interest is the distributions of ther- 
mal and CR pressures inside groups/clusters. The differ- 
ence in the distributions of the two pressure components 
is an important detail, because, according to the condition 
for hydrostatic equilibrium, 

dP tot {r) _ GM cl (r)p gas {r) ^ 
dr r 2 
so, it is the total pressure gradient that responds to the 
group/cluster mass enclosed in a volume of radius r. Thus 
the spatial distribution (gradient) of P cr is as important 
as the amount of CR energy content (E cr ~ 3 P cr ) itself, 
once P cr becomes dynamically significant. 

First we consider the ratio of CR to thermal pres- 
sure evaluated within each computational cell inside 
groups/clusters; i.e., (P cr / Pth) ceil ■ The average of this 
quantity over the cells inside each group/cluster volume 
is plotted in the left panel of Fig. 2, as a function of the 
group/cluster temperature. This is similar, but not iden- 
tical to (P cr / Pth)d in Fig. 1, The standard deviation of 
(P cr / Pth) ceil values within each group/cluster is shown in 
similar fashion on the right panel of the same figure. It 
is clear that the dispersion around the average value is as 
large as the average itself, indicating a strong variation of 
(Per I Pth) ceil inside each group/cluster volume. Note that 
the gas temperature and the slope of the CR distribution, 
q, are approximately uniform within groups/clusters. So, 
this pressure behavior should be a reflection of the differ- 
ent spatial distributions of gas and CRs in the simulation. 

In order to quantify the difference in spatial distribu- 
tions of the two pressures, we define a pressure- weighted 
mean square radius 

R] = (3-5) 

where Pi and rj indicate the pressure (either thermal or 
CR) and the distance from the group/cluster center of the 
i-th computational cell, respectively. In Fig. 3 we plot 
the ratio of Ri relative to the CR and thermal pressure, 
i.e., Ri(p cr )/Rj( Pth y The plot shows that this ratio is close 
to one, with a marked tendency to values slightly larger 
than one. That is, CR pressure would be distributed more 
diffusely than gas pressure in groups/clusters. Caution is 
needed here, since the diameter of the collapsed objects 
covers only about 5 computational cells. Thus, although 
the collapsed objects have been formed with adequate res- 
olution to assure their basic properties, the fine details 
of their structures have not been captured. Nevertheless, 
the systematic difference could be connected to the mech- 
anism of CR production. In fact, strong shocks in the 
simulations arc more commonly located at the outskirts 
of the collapsed object. There the ratio of CR to thermal 
pressure is therefore higher, causing i?/( Pcr ) to be slightly 
larger than Ru Pth y 

3.3. Gamma Ray Flux 

A direct observational consequence of CR protons in 
groups/clusters is the 7-ray emission from n° decay. The 
7-ray emissivity was calculated in each cell from the gas 
density and the CR proton distribution. The cross sec- 
tions were computed according to the GALPROP rou- 
tines (Moskalenko & Strong 1998). The number of n° pro- 
duced in each hadronic interaction, £ n o, increases rather 
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slowly with the proton energy, E p , roughly as ~ 
[(E p - S t/ls )/GeV]V4 for E ths < E p < 10 4 GeV (E ths = 
1.22 GeV is the energy threshold of the process; see, 
e.g., Mannheim & Schlickeiser 1994). For a proton power- 
law distribution with kinetic energy, f C r(T p ) oc T~ q+2 , 
CRs at energy T p generate a number of ir° which roughly 
scales as oc (Tp/GeV)-^ 5 / 2 ) . Then, since q > 4 the 
majority of the integrated 7-ray flux is contributed by CR 
protons in the low energy component. 

In Fig. 4 we report the expected 7-ray flux above 100 
MeV, Fj, from a volume within 1.3 /i Mpc from the 
group/cluster center, as a function of the core temper- 
ature, T x . An integration volume larger than a typical 
group/cluster core region of 0.5/i _1 Mpc is chosen, because 
the CR proton distribution extends out further to where 
the accretion shocks are found. When fitting the F 7 — T x 
relation with a power law curve from a simple \ 2 analysis 
we get 

, rp \ 2.95 

F 1 = 7 A x 10~ 9 ( - 72 ^ v J counts s" 1 cm" 2 . (3-6) 

We note here that the spread about the average at a given 
T x exhibited in Fig. 4 is of order of a few. This scatter, also 
pointed out in §3.1, is almost certainly real and is a reflec- 
tion of the different peculiar formation history and current 
dynamical state that can characterize a group/cluster with 
a given temperature. With this relation we find that for 
a Coma-like cluster with temperature of T x — 8.3 keV the 
mean 7-ray flux would be <~ 1.4 x 10~ 8 counts s _1 cm~ 2 . 
Similarly, after rescaling the flux for the appropriate dis- 
tance, we can compute F 7 expected for other clusters of 
known temperature. Thus, we find F 7 ~ 9.8 x 10~ 9 for 
a temperature of 6.3 keV and a distance of 55 /i _1 Mpc 
characteristic of the Perseus cluster (Schwarz et al. 1992); 
and Fj <~ 3.8 x 10~ 9 for a temperature of 1.8 keV and 
a distance of 14 /i _1 Mpc as estimated for the M87- Virgo 
system (Bohringer et al. 1994). 

The values found above for nearby clusters are well be- 
low the limit set by the EGRET experiment of 4 x 10~ 8 
counts s _1 (Sreekumar et al. 1996). However, they should 
be easily detectable by GLAST, with a sensitivity an order 
of magnitude below the above value. Our estimate of the 
7-ray fluxes is somewhat lower than the values computed 
by other authors (Dar & Shaviv 1995; EnBlin et al. 1997), 
which are close to, or slightly in excess of, the EGRET 
upper limits. This comparison is so even after correcting 
our results for the aforementioned systematic underesti- 
mate due to resolution effects. Differences in estimates of 
the 7-ray flux are expected, given the numerous differences 
in our physical assumptions and methodology when com- 
pared to the previous authors. For example, the spatial 
distribution of CRs given by the simulation in our case 
was assumed, by contrast, to be uniform in Dar & Shaviv 
(1995), or such as to produce a constant ratio of CR to 
thermal pressure in EnBlin et al. (1997). In addition, we 
have used a fixed emitting volume within a radius of 1.3 
/i _1 Mpc, whereas the previous authors used ^4 /i _1 Mpc 
for the size of the Coma cluster (Dar & Shaviv 1995; EnBlin 
et al. 1997), and ^1.3 /i _1 Mpc for a Virgo-like type of 
cluster (Dar & Shaviv 1995). Also, for the slope of CR en- 
ergy distribution those authors (Dar & Shaviv 1995; EnBlin 
et al. 1997) borrowed the empirical value from the Galac- 



tic case; i.e., assumed q = 4.7, unlike q = 4.0 — 4.2 from 
our simulation. 

Our estimate of the 7-ray flux for the Coma cluster is, 
on the other hand, compatible with the value computed by 
Colafrancesco & Blasi (1998), although our temperature 
dependence in Eq. (3.7) is quite a bit steeper than the one 
they presented. They computed the emission within the 
virial radius, whereas a fixed volume has been used in our 
estimate. In addition, they accounted analytically for a 
weak phenomcnological dependence of the baryon fraction 
on the cluster size, while the baryon density from numer- 
ical simulation has been used in our estimate. We note 
that the expected functional form of the 7-ray flux can be 
modeled as 

F 1 (xN cr n b , (3-7) 

where N cr is the total number of CR protons and nj, is the 
average group/cluster baryonic mass density, both inside 
a fixed radius. So, most of the temperature scaling in our 
estimate is accounted for by the following facts: (1) The 
kinetic energy power available for CRs is proportional to 
T 2 (see Miniati et al. 2000). That is, E cr oc T 2 . On the 
other hand, with a constant momentum slope in the CR 
distribution function, N cr oc E cr . Hence, N cr oc T x holds 
approximately in the simulation. (2) The mean mass den- 
sity inside a fixed volume (thus, also the associated baryon 
mass) scales almost linearly with the temperature, which is 
compatible with observations (e.g., Edge & Steward 1991; 
Mohr et al. 1999). Together, those explain the origin of 
the ~ X! 3 dependence found in the simulation presented 
here. 

One important point of our findings is that the calcu- 
lated values for the 7-ray flux are well below the upper 
limits set by EGRET, even while a large fraction of the 
total energy in the ICM gas is, in fact, stored in CRs. 
This result differs from Blasi (1999), who finds that the 
amount of energy in the CR component must be well be- 
low the equipartition value in order not to violate the same 
observational limits (except for the extreme case of a uni- 
form distribution of CRs throughout clusters). Blasi's re- 
sult derives from his adoption of a central point source for 
cluster CRs which then must diffuse throughout the clus- 
ter. That leads to a CR distribution more concentrated 
towards higher thermal gas densities in the core than if 
the CRs are produced by structure shocks as in our case. 
Since in a denser environment the CRs experience many 
more interactions, a higher 7-ray flux is expected. 

Finally, we have found that there is a tight correlation 
between the CR pressure and the 7-ray flux that can be 
fit by 

f p xO.64 

P cr = 2.7X10- 11 — — erg cm" 3 . 

\ 10 M counts s cm z / 

(3-8) 

This is shown in Fig. 5. The above scaling is compati- 
ble with F 7 oc T 3 and P cr oc T 2 , if the slope of the CR 
momentum distribution varies only slightly (as it is the 
case here). P cr — F 7 is manifestly a cleaner relation than 
P cr — T x (not shown in this paper). The outcome is not 
surprising and is due to the direct physical relationship be- 
tween P cr and Fry, i.e., F 7 oc j n gas P cr (for fixed q). Thus 
the P cr — F 7 relation probably allows the best determina- 
tion of the amount of CR pressure in groups/clusters. In 
addition, measuring the F 7 — T x relation and comparing 
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it to the numerical results will provide an important test 
for our numerical treatment and general understanding of 
CR injection, transport, and acceleration in group/cluster 
environment. 

3.4. X-ray and "{-ray Images 

In Fig. 6, we present synthetic images of 7— ray emis- 
sion from 7T° decay (left) and X-ray emission from ther- 
mal bremsstrahlung (right) for a rich group of galaxies 
in the simulation with T x ~ 3 keV. A larger sample of 
groups/clusters from the numerical simulation was stud- 
ied by Miniati (2000), but their main properties can also 
be summarized by the findings below. The volume consid- 
ered for the realization of the images in Fig. 6 is a cube 
about 14 h~ 1 Mpc on a side, centered on the collapsed ob- 
ject. Each synthetic map has a grey-scale bar indicating 
the logarithmic value of the imaged quantity. The physi- 
cal units are in 'erg cm -2 per pixel' for the flux density of 
X-rays, and in 'counts cm" 2 per pixel' for the flux density 
of 7— rays (see §2.3 for more details). In addition to the 
synthetic images, in Fig. 7 we also present, for the same 
object, two-dimensional slices of the following quantities: 
the CR proton number density in units of cm -3 (top left), 
the velocity field (top right), the contours of shock com- 
pression (V • v - bottom left) and the gas number density 
in units of cm -3 (bottom right). The slices are through 
the object center and perpendicular to the line of sight of 
the synthetic images. The side of the images in each panel 
of Fig. 7 is about 5.3 /i _1 Mpc. 

First, we note that the synthetic 7-ray image exhibits an 
irregular morphology, somewhat different from the smooth 
slightly ellipsoidal shape of the thermal X-ray image. This 
is due to the fact that the CRs, responsible for the 7-ray 
emissivity, are sensitive to the particular shock distribu- 
tion in the ICM. When even a mild shock crosses through 
a group/cluster, the injection of fresh particles over a rel- 
atively short time can significantly enhance the CR pop- 
ulation. This should definitely affect the shape of the 7- 
ray image. But it is not so for the X-ray image, because 
the effect of a weak shock is only a modest increase of 
the density (square) on which the emissivity primarily de- 
pends. That effect is likely to be blended away after line of 
sight integration (although finite resolution may limit the 
amount of visible features as well). In fact, closer inspec- 
tion by means of two-dimensional slice images shows that 
enhanced cosmic ray density occurs downstream of "inter- 
nal shocks"; i.e., post-shock flows in the central regions 
of groups/clusters (Miniati 2000). In Fig. 7 the verti- 
cally elongated, high CR density structure (top-left panel) 
is enclosed by a Mach surface where a supersonic flow is 
suddenly decelerated by a shock. This region corresponds 
to a vertical dent in the velocity vector field (top-right) in 
the direction N-W from the center of the panel. Note that 
the cosmic rays are where most of the intra-cluster gas is 
located (bottom-right panel) in addition being near the 
shock, and, therefore, where the injection rate is higher. 
Such higher level of structure and more irregular distribu- 
tion of the CRs, as compared to the gas density, explains 
in part the high pixel to pixel fluctuation of the ratio of 
CR to thermal pressure that was found in §3.2. 

Finally we note that the regions of low surface bright- 
ness corresponding to the same factor below the peak value 



are slightly more extended in the 7-ray than in the X-ray 
image in accord with findings in §3.2. Also, the highest 
7-ray surface brightness appears more concentrated than 
the X-ray brightness distribution. The X-ray emissivity 
is proportional to the square of the gas density, whereas 
the 7-ray emissivity is proportional to the product of the 
gas density with the CR density. Thus, this result sim- 
ply means that the CR protons are slightly more concen- 
trated by number than the gas in the group/cluster core 
region. This result does not contradict the finding in 3.2 
that the CR pressure distribution is less concentrated than 
gas pressure. Rather, it indicates that adiabatic compres- 
sion is effective and has reduced the ratio (P cr /Pth) m the 
center of the collapsed objects. 

4. DISCUSSION 

As we have shown in the previous section, a significant 
fraction of the total energy associated with baryons in- 
side a group/cluster could be stored in cosmic rays as a 
consequence of diffusive particle acceleration at structure 
formation shocks. This fact bears important consequences 
that we will try to address in the following discussion. 

Firstly, if the pressure provided by cosmic rays, P cr = 
E cr /3, is large enough, it can affect the dynamics of 
the intra-cluster medium and, therefore, both its evolu- 
tion and equilibrium. This is of great concern because 
groups/clusters of galaxies are invaluable probes to test 
cosmological theories and to measure key cosmological pa- 
rameters (e.g., Bahcall et al. 1999, and references therein). 
In fact the present-day abundance (number density) of rich 
clusters of galaxies sets a strong constraint on the total 
mass content of the universe and the normalization of the 
power spectrum of the density perturbation by imposing 

a 8 nU 2 ^ 0.5 ± 0.05 (Bahcall & Cen 1992; White et al. 
1993a; Eke et al. 1996; Viana & Liddle 1996; Pen 1998). 
In addition, the evolution of the number density of rich 
clusters allows one to break the degeneracy of the above 
result and is used to determine both as and O m (Carlberg 
et al. 1997; Bahcall et al. 1997). It is not clear whether and 
how the evolution of structure would be affected by a non- 
thermal dynamical component. Clearly, since most of the 
mass is dark, the growth of the density perturbation would 
be unchanged for the most part. However, since the ob- 
servable universe is made of baryons, the specific processes 
that determine their dynamical and thermal evolution are 
of crucial importance, as they also greatly affect cluster 
observational properties. In this respect the effects pro- 
duced by the cosmic ray component could be important, 
even though the underlying large scale structure remains 
unaffected. 

Furthermore, a substantial cosmic ray pressure compo- 
nent could contribute to the dynamical support of the 
intra-cluster medium against gravitational collapse. This 
would affect the estimate of the total cluster mass derived 
from observations and, in turn, both of the baryonic frac- 
tion there and of the total mass of the universe (White 
et al. 1993b). Results from a number of studies have 
suggested that mass estimates based on the hydrostatic 
equilibrium assumption and X-ray measurements tend to 
be somewhat smaller than those derived from virial esti- 
mates and gravitational lensing (Markevitch & Vikhlinin 
1997; Horner et al. 1999; Nevalainen et al. 2000; Roussel 
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et al. 2000; Miralda-Escude & Babul 1995; Wu 2000, see 
also Miniati 2000 for an extensive discussion on the issue). 
Part of this could well be the consequence of dynamical 
effects due to the cosmic-ray pressure and also the mag- 
netic field pressure. However, the mass discrepancy issue 
is still controversial and the precision of the current mea- 
surements, at the level of 20% accuracy, does not allow a 
strongly conclusive statement at this point. As discussed 
in the previous section, 7-ray observations appear promis- 
ing in this respect, since according to our prediction the 
expected 7-ray flux from Coma-like clusters should be well 
above the detection threshold of GLAST (see also Blasi 
1999; Dolag & Enfilin 2000). 

From the theoretical side, recent numerical simulations 
have shown that, in order to construct a viable and re- 
alistic depiction of the ICM, the various processes taking 
place there need to be accounted for in sufficient detail. In 
particular, it has become clear that the effect of radiative 
cooling in cluster cores produces significant quantitative 
differences in numerical simulations in which it is allowed 
(e.g., Katz & White 1993; Suginohara & Ostriker 1998; 
Pearce et al. 2000; Lewis et al. 2000). In particular Sugi- 
nohara & Ostriker (1998) found that cooling can become 
catastrophic in the cluster cores unless prevented by some 
additional physical processes. Similarly Lewis et al. (2000) 
found that radiative cooling can have dynamical effects on 
the cluster structure and evolution. In particular, they 
concluded that the consequences of cooling are global and 
affect the cluster as a whole, despite the fact that strong 
cooling is localized in the central region of a cluster. Ac- 
cording to the study by Lewis et al. (2000), however, the 
catastrophic character emerging in the Sughinohara & Os- 
triker's simulations is largely inhibited by the feedback of 
star- formation (gas removal and heating). Nevertheless, 
that does not solve the cooling problem completely. In 
fact, the star formation ensuing from the cooling of the 
gas produces too large a stellar component (30% of all the 
baryons instead of the observed 10% fraction), too high an 
X-ray luminosity by a factor ~ 3 and too high a velocity 
dispersion (Lewis et al. 2000). These excesses are driven 
by a very high density, stellar dominated core resulting 
from the effect of cooling. In this respect, the presence of 
a significant CR component could reduce the overly dra- 
matic effect of radiative cooling and recover some of the 
observed cluster properties, at least in two ways. Firstly, 
CRs provide an additional non-thermal pressure, which is 
not dissipated by radiative effects. This hinders the con- 
traction of the cooling gas; therefore, prolonging the cool- 
ing time and decreasing the rate of conversion of gas into 
stars. Secondly, low energy CR ions provide a source of 
heating that tends to balance cooling, once again soften- 
ing the effect of the latter (Rephaeli 1977). In this respect, 
CRs are also likely to affect the dynamics of a cooling flow 
by means of the two generic mechanisms described above. 

5. CONCLUSIONS 

We have carried out a computational study of produc- 
tion of CR protons at cosmological shocks associated with 
the large scale structure in a SCDM universe. We have 
achieved this by carrying out the first numerical simulation 
of structure formation that includes directly shock acceler- 
ation (in the test-particle limit approximation), transport 



and energy losses of the CRs. CR injection takes place 
at shocks according to the thermal leakage prescription, 
leading to the injection as CR of a fraction about 10~ 4 of 
the thermal protons passing through a shock. According 
to our results, cosmic ray ions may provide a significant 
fraction of the total pressure in the intra-cluster medium. 
The conclusion cannot be made strictly quantitative yet, 
because the complex physics regulating the acceleration 
mechanism cannot be fully simulated, and our simulated 
group/cluster structures are still rather coarse. However, 
we expect the CR pressure may account for a few tens 
of percents of the total ICM pressure. The cosmological 
consequences of this result were addressed in the previous 
section. 

A major step forward will be made possible by the ad- 
vent of the next generation of 7-ray facilities, i.e., GLAST. 
In fact, we expect 7-rays will be detected for relatively 
nearby massive clusters. That development will probe 
directly the cosmic ray content in clusters of galaxies 
(see §3.3). In addition, 7-ray imaging and spectroscopy 
will enable us to infer the spatial distribution of the CR 
density and pressure, once the gas distribution is known 
(e.g., through X-ray data). That will translate in direct 
information on the the nature of the CR sources. In fact, 
if most of the CRs have been expelled by active galax- 
ies, then their distribution would not be as widespread as 
in the case where the primary sources are cosmic shocks. 
This adds to the wealth of critical information provided by 
observation in this band. However, most probably only the 
largest clusters and only their innermost regions of high- 
est emission will be probed by these instruments, owing to 
the very low surface brightness in the 7-ray band. Never- 
theless, those detections would still be invaluable for the 
study and a much deeper understanding of the dynamics 
of these objects within the next few years. 
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Fig. 1. — Ratio of CR to thermal pressure averaged over the group/cluster volume within 0.5 h 1 Mpc plotted as a function of group/cluster 
core temperature. 
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Fig. 2. — Left panel: average of the cell by cell evaluation of (Per / Pth) cell ■ Right panel: standard deviation. 
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Fig. 3. — Ratio of pressure-weighted rms radii of CR to thermal pressure defined in Eq. (3.5). 
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Fig. 4.— 



7-ray flux as a function of group/cluster core temperature. 
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Fig. 5. — CR pressure as a function of 7-ray flux. 
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Fig. 6. — Synthetic images in 7-rays from ir° decay in units counts s 1 cm 2 per pixel (left) and X-ray from thermal bremsstrahlung in 
units erg s~ 1 cur 2 per pixel (right) from a cosmic volume of (14 ft _1 Mpc) 3 . 
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Fig. 7. — Two-dimensional slice maps of CR proton number density in units cm 3 (top left), velocity field on that plane (top right), 
contours of shock compression (V ■ v - bottom left) and gas number density, again in units cm~ 3 (bottom right). The side of images is 5.3 
ft-iMpc. 
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